Mixing property of triangular billiards 
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We present numerical evidence which strongly suggests 
that irrational triangular billiards (all angles irrational with 
7r) are mixing. Since these systems are known to have zero 
Kolmogorov-Sinai entropy, they may play an important role 
in understanding the statistical relaxation process. 

PACS number: 05.45.+b 



After the pioneering paper of Fermi, Pasta and Ulam 
jjj and the mathematical works of Kolmogorov school, 
the modern ergodic theory can now account for the rich 
variety of different statistical behaviours of classical dy- 
namical systems. These properties range from complete 
integrability to deterministic chaos. The relatively few 
rigorous results so far available have provided a firm 
guide for a large amount of analytical and numerical work 
which gave a strong impulse to the field of nonlinear dy- 
namics. Needless to say, several problems remain to be 
solved. For example it is known that mixing property 
guarantees correlations decay and relaxation to statistical 
equilibrium. It is however not known whether or not this 
property is sufficient for a meaningful statistical descrip- 
tion of the relaxation process or whether the strongest 
property of positive KS entropy is required. Even much 
less clear is the situation in quantum mechanics. For 
example, to what extent the statistical distribution of 
energy levels of a quantum conservative hamiltonian sys- 
tem, is related to the different properties in the ergodic 
hierarchy is, as far as we know, an open question. 

Examples of classically completely integrable or deter- 
ministic random systems have been widely studied in the 
literature. However, to our knowledge, there are no phys- 
ical examples of systems which possess the mixing prop- 
erty only (with zero KS entropy). In this respect, the best 
candidates are billiards in 2-d triangles but, in spite of 30 
years of investigations, no definite statement can be made 
concerning their dynamical properties |Q. Moreover, || 
"a prevailing opinion in the mathematical community is 
that polygonal billiards are never mixing, but this has 
not been established" . 

In this paper we provide strong numerical evidence 
that generic triangular billiards, with all angles irrational 
with it, are mixing. This result can play an important 
role in the understanding of non-equilibrium statistical 
mechanics since it fills a gap in the ergodic hierarchy. 
Moreover, since the local dynamical instability in these 
systems is only linear, they have zero algorithmic com- 



plexity. This means that, even though it may be very dif- 
ficult, their analytical solution is not impossible in prin- 
ciple and this may prove to be very important for the 
future development of nonlinear dynamics. 

We consider the motion of a point particle, with unit 
velocity, inside a generic triangular billiard (with all an- 
gles irrational with tt, in general). Let us introduce carte- 
sian coordinates (x,y), and align one side of the triangle 
along the x— axis (y = 0), starting at the origin x = 0, 
and let us choose its length to be unity. We denote the 
angle at the origin (0,0) by a and the angle at (1,0) by 
/3. Hence the other two sides of the triangle are given by 
the eqs. y = (tana)cc, and y = (tan/3)(l — x). 

We would like to draw the attention to the following 
fact: the dynamics of a point particle in a triangular 
billiard is equivalent to the motion of three particles on 
a ring with different masses, mi, mi, and 7713. Let the 
ring have circumference 1. The motion of 3-particles Id 
gas with coordinates qi,q 2 ,q 3 (qi < 32 < 33 < Qi + 1) 
is governed by the hamiltonian H = ^miq\ + ^miq\ + 
\m 3 q^, with elastic collisions at qi = qi, qi = q 3 , and 
93 = Qi + 1- Introducing the notation 

f = [y/mlqi , y/rn^qi , y /m 3 'q 3 ) 

and the orthogonal transformation 
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with M = mi + m2 + ^3, the hamiltonian writes: 
H = i(i 2 +y 2 + i 2 ) where, as can be checked by straight- 
forward calculations, the motion in the (x, y)— plane is 
bounded by specular reflections from the sides of the tri- 
angle with angles 



tana 



m 2 M 
mim 3 ' 



tan/3 



miM 
m 3 m 2 ' 



tan 7 



m 3 M 
mimi 

(1) 



The motion in the z~ coordinate trivially separates and 
corresponds to the center-of-mass motion in the 3- 
particles Id gas. 

Notice, however, that triangular billiards with one an- 
gle larger than it/ 2 cannot be related to the 3-particles 
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Id gas, since from the above formula it follows that 
a, (3, 7 < 7r/2. A particular case is given by the isosceles 
triangle which is dynamically equivalent to the right tri- 
angle billiard, with e.g. 7 = tt/2. The latter corresponds 
to 7713 = 00, i.e. to the Id motion of two particles with 
masses mi , 777,2 between hard walls Q . 

The dynamics of triangular billiards can be divided 
into three classes: 

(A) All angles rational with it. The dynamics of such 
triangles is not ergodic; in fact it is pseudo-integrable, 
i.e. it possesses 2d invariant surfaces of high genus in 4d 
phase space. 

(B) Only one angle rational with tt. Such are generic right 
triangles, for example. In recent numerical experiments 
evidence has been given that right, irrational, triangular 
billiards are ergodic and weakly mixing [l] . 

(C) All angles irrational with tt. Surprisingly, to the best 
of our knowledge, this generic class of triangles has been 
somehow overlooked in previous numerical studies and 
will be the main object of the present paper. It is within 
this class that one may now hopefully look for ergodic 
and mixing behaviour. 

From the rigorous point of view not much is known be- 
yond the fact that the set of ergodic triangles is a dense 
G$ (intersection of a countable number of dense open 
sets) in a suitable topology. We recall that a dynami- 
cal system T : x(0) — > x(t) with invariant measure n(x) 
in phase space A4 is mixing if, for any L 2 pair of ob- 
servables in phase space, their time correlation function 
asymptotically vanishes 
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The map T may represent a continuous flow (t real) or 
a discrete map (t integer). In this paper we shall mainly 
consider the dynamics given by a discrete Poincare map 
which corresponds to the collisions of the orbits with the 
horizontal side y — 0. The reduced phase space — surface 
of section (SOS) — is a rectangle, parametrized by the 
coordinate < x < 1 and by the corresponding canonical 
momentum — 1 < p x — sin$ < 1 (# = angle of incidence), 
with the invariant measure d\i = dxdp x . 

As a first step we have performed an accurate check of 
ergodicity. We have done this in two independent ways: 

(i) By dividing the phase space (SOS) in a large num- 
ber N = N\ x N\ of cells and then computing the number 
n(t) of cells which are visited by a single orbit up to dis- 
crete time t I9I. Then we computed the average relative 
measure r(t) = (n(t)/N) of visited phase space up to 
time t where (.) denotes phase space average over suffi- 
ciently many randomly chosen initial conditions. As it 
is known, for the so-called random model of completely 
uniform (ergodic) and random dynamics, one can derive 
the simple scaling law [p| 



Expression (|^) (for arbitrary but sufficiently fine mesh 
N) should be considered as a sufficient but not necessary 
condition for ergodicity. In fact, the dynamics would 
obey the law (|2|) only when the system does not possess 
any non-trivial time scale. In Fig.l we show the results 
of numerical computations for triangles of class (B) and 
(C). We take for example a right triangle B with angles 
a = (\fb — l)7r/4, (3 — tt/2 — a, 7 = tt/2 and a generic 
triangle C with angles a — (\/2 — 1)tt/2 : (3 = {-\fb — 
1)71-/4,7 = 7T — a — /3. It is seen that there is a drastic 
difference between the two cases: the generic triangle C 
excellently follows the law (^), whereas the right triangle 
B, even though ergodic, Q strongly deviates and explores 
the phase space extremely slowly, with a hierarchy of 
long time-scales related to a strong sticking of the orbit 
in momentum {p x ) space. 




FIG. 1. The relative measure r(t) of the visited phase 
space as a function of discrete time t. The full and dashed 
curves refer to different discretizations of phase space, namely 
to N — 10 6 and N = 10 4 cells, respectively. The two upper 
curves refer to the generic triangle C whereas the two lower 
curves refer to the golden right triangle B (see text for details) . 
The random model, expr.(2), is given by the dotted-dashed 
curve which is almost indistinguishable from the numerical 
curve for the triangle C at N = 10 6 . In order to better display 
the excellent agreement with expr.(2), we plot in the inset 
the function 1 — r(t) for the generic triangle C in the semi 
log-scale. At N = 10 6 , and N = 10 4 , averages over 200, 
and 4000 orbits with randomized initial conditions have been 
used, respectively. 

(ii) By comparing the time averaged correlation func- 
tion for a single, but very long orbit, C*(i) = 

limT-*oo h T^n = \ a(n)a(n + t) with the phase averaged 
correlation function C%(t) = (a(0)a(t)) which is com- 
puted by monte-carlo averaging over many short orbits 
with different randomized initial conditions. In Fig. 2 we 
plot the difference between the two curves for the generic 
triangle C which turns out to be of the same order as the 
statistical error. 



r(i) = r RM (t) = 1 - expH/AO- 
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FIG. 2. The velocity autocorrelation function Cp x (t) in 
log- log scale for the generic triangle C (full curve), which 
is obtained as the average over 1.2 • 10 6 different orbits of 
length 32768. The dashed line gives the slope -1.00. The 
dotted curve gives the estimated statistical error of the cor- 
relation function, while the dots mark the difference between 
the phase averaged and the time averaged correlation func- 
tion \Cp (t) — Cp x (t)\ where the time average is computed 
from a single orbit of length T = 32768 ■ 1.2 • 10 6 with ini- 
tial condition xo = 0.23456, p x o = 0.34567. The statistical 
error is estimated as the standard deviation of a sequence of 
M = 1000 partial averages of length T/M each. 

The above numerical results are surprisingly much more 
clear than expected; they demonstrate ergodic behaviour 
and they make reasonable to expect mixing behaviour 
also. We now turn our attention to the latter question. 

We have performed extensive numerical computa- 
tions of autocorrelation functions of different observables, 
namely, momentum (horizontal component of velocity) 
v = p x , symmetrized position x 1 = 2x — 1, and also 
characteristic functions of sets in phase space. In all 
cases we have found clear numerical evidence of power- 
law decay of correlation functions over about 4 orders of 
magnitude. For a given triangle, correlation functions 
of different observables, decay with the same empirical 
exponent which is typically very close to — 1. In Fig. 2 
we show the velocity correlation function of the trian- 
gle C, which decays with exponent with the fitted value 
— a = —0.90 ± 0.02. However, since the triangle C is 
not too far from the right triangle we have chosen an- 
other, similar generic triangle D with the irrational an- 
gles a = (s/2 — 1)tt/2, /3=1,7 = 7t — a — /3 and show in 
Fig. 3 its velocity and position correlation function which 
decay with empirical exponent — a = —0.94 ±0.04. We 
have performed extensive numerical calculations of corre- 
lation functions for many different generic triangles, and 
the exponent of decay have always been very close to 
— 1. There are even cases where the fitted decay expo- 
nent seems to be slightly below —1 (a > 1), however in 
such cases, correlation functions appear much more noisy 
and it is difficult to make precise statements. 
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FIG. 3. The velocity (a) and position (b) correla- 
tion function, Cp(t) and C*, (t), for an orbit of length 
10 11 in the generic triangular billiard D. Initial condition: 
x = 0.23456, pso = 0.34567. 



Finally, we have performed a different and powerful 
statistical test, namely we have computed the Poincare' 
recurrences or the return time statistics, i.e. the proba- 
bility P(t) for an orbit not to stay outside a given subset 
A C M. for a time longer than t. It has been conjec- 
tured HQ (see also || for accurate definition and fur- 
ther references) that the integrated probability Pi(t) = 
J2t?=t should be intimately connected to the corre- 
lation decay. More precisely, if the correlation function 
decays as a power-law with exponent —a, then the in- 
tegrated return probability should decay with a similar 
exponent P(t) ~ i - " -1 , P„(i) ~ t~^, where fi = a. In 
Fig. 4 we show the recurrence probability P(t) and the 
integrated probability Pi(t) computed w.r.t. half-space 
A = {(x,p x );Px > 0} for the same data as in Figs. 2, 3 
(triangular billiard C). We have found that numerical re- 
sults are perfectly consistent with fj, = a = 1, i.e. with 
1/t 2 decay of return probability P(t). 

To summarize, the numerical obtained value of expo- 
nent a is close to -1 but in some cases is slightly less 
and in other cases slightly larger. We attribute this fact 
to long transient times which are caused by the intricate 
interplay of arithmetic properties of angles a, /3, 7 and 
by the presence of periodic orbits. In this connection we 
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would like to mention that we have numerically investi- 
gated the existence of periodic orbits and we have found 
that their number increases very slowly with their pe- 
riod. The number of non-equivalent periodic orbits with 
lengths up to I (collisions with the boundary) typically 
increases slower than I. 
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FIG. 4. The return probability P(t) (full curve), and the 
integrated return probability Pi(t) (dashed curve), for the 
same orbit (triangle D) as in Fig. 3. The dashed-dotted lines 
have slopes —2 and —1 respectively. 

A qualitative argument which leads to the 1/t decay 
of correlations reads as follows: let us consider a fixed, 
sufficiently small, region A in phase space around a peri- 
odic orbit and study the I iterates of the Poincare map, 
where I is the period of the orbit. Due to the linear 
(in)stability, the orbit which starts in a small square a e 
of side e centered on the periodic orbit and contained in 
A, will remain in A for a time t > 1/e. It follows that 
the probability P(l/e) for an orbit to stay in the region 
A for a time t > 1/e is proportional to the probability 
of the orbit to visit the square a e . Now, the linear in- 
stability implies that an orbit which enters the square a e 
will remain inside a e for a constant time which does not 
depend on e as e — > 0. As a consequence, from the Liou- 
ville theorem it follows that the probability for an orbit 
to visit the square a e will be proportional to its area e 2 . 
This is the probability P(l/e) to remain in the region A 
for a time t> 1/e. This leads to the relation P{t) oc 1/t 2 
as confirmed by numerical results. 

Since the question discussed in this paper is a very del- 
icate one, we have put particular attention to the accu- 
racy of our numerical computations. Among the several 
different tests, we have developed two independent com- 
puter codes, one based on the billiards dynamics in a 2d 
plane, and the other based on the three particles Id gas 
dynamics: the results agree. Moreover, since the instabil- 
ity here is only linear in time, and the machine precision 
is ~ 10~ 16 , the numerical errors are always below the 
statistical errors. 

The central question is to what extent the results pre- 
sented here can be considered as a definite or convincing 



evidence for the mixing property of generic triangular bil- 
liards. The correlation decay of some particular class of 
functions is certainly compatible with a weaker property 
than mixing. However, in this paper we have compared 
the correlations decay of different variables (velocity, po- 
sition, characteristic functions of phase space sets etc). 
We have also considered the dynamics given by the dis- 
crete Poincare map relative to different sides of the same 
triangle: all these correlations exhibit a power law decay 
with exponent very close to —I, and make very plausi- 
ble the conjecture that correlation functions in generic 
triangles decay as oc 1/t. Of particular significance is 
the behaviour of the return probability P(t) which de- 
cays with the expected power law cx 1/t 2 . All the above 
results (including those described in Fig.l) lead to the 
conclusion that, outside any reasonable doubt, generic 
irrational triangles are mixing. The particular case of 
isosceles or right triangles is much less clear. The be- 
haviour of correlations in such case is noisy and the ex- 
ponent of the decay is too small to allow for any definite 
conclusion. 

The fact that generic triangles have zero Kolmogorov- 
Sinai entropy and yet they have a very nice mixing be- 
haviour without any time scale, may prove to be very 
useful for understanding the dynamical basis of the re- 
laxation process to the statistical equilibrium. Moreover, 
the analysis of their quantum behaviour will contribute 
to the current efforts for the construction of a statistical 
theory of quantum dynamical systems. 
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